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Abstract 

Many natural populations are well modelled through time-inhomogeneous stochastic pro- 
cesses. Such processes have been analysed in the physical sciences using a method based on 
Lie algebras, but this methodology is not widely used for models with ecological, medical and 
social applications. This paper presents the Lie algebraic method, and applies it to three 
biologically well motivated examples. The result of this is a solution form that is often highly 
computationally advantageous. 

1 Introduction 

Stochastic models based on Markov chains are important in many ecological, medical and social 
contexts. In these contexts, where populations are modelled, there are often external influences 
that act on the system in a manner that varies over time, leading to a time-inhomogeneous Markov 
chain model and corresponding technical difficulties for analysis [10] , 

Wei and Norman [14] proposed a method for dealing analytically with time-inhomogeneous 
Markov chains, based on Lie algebraic methods. The idea of combining Lie algebras and symmetry 
considerations with Markov chains has continued to attract theoretical interest in a variety of 
contexts QUEUSE!!- 

At the same time, there is a more applied desire to have numerically efficient methods to analyse 
Markov-chain population models, one option for which is the use of matrix exponentials [8j HOj . The 
aim of this paper is to explain how Lie algebraic methods can be used to derive matrix exponential 
solutions to time-inhomogeneous Markov chains that are applicable to population modelling. In 
contrast to other applications, symmetries of these systems are not a guide to the appropriate Lie 
algebra to use in solution of population models; a certain amount of trial and error is necessary. 
The focus of this paper is therefore on three examples in population modelling where it is possible 
to define an appropriate Lie algebra, and a discussion of the potential benefits of doing so. 

2 Methodology 
2.1 Lie algebras 

In general, a Lie algebra over a field F is an F-vector space V, together with a bilinear map called 
the Lie bracket. Elements of the vector space are written u,v,w S V . The Lie bracket is written 
[it, v] G V and obeys 

[u,u} = 0, [u, [v, w]} + [v, [w, it]] + [w, [u, v]} = . (1) 

We will be interested in the vector space GL(n, R), i.e. the set of real- valued n x n matrices, and 
will define the Lie bracket for A, Y g GL(n, R) through the commutator 

[A, Y] := XY — YX , (2) 

which can be readily seen to satisfy ([T]). It is often convenient to define an adjoint endomorphism 
operator, ad, to represent the Lie bracket: 

(adA)F :=[A,F] , (3) 



f 



so that multiple applications of the Lie bracket can be concisely written as e.g. (a.dX) 2 Y = 
[X,[X,Y]] . 

2.2 Time-inhomogeneous Markov processes 

Suppose that pit) is a probability vector, i.e. a vector with values p n (t) > summing to unity 
that represent the probability that an integer stochastic variable takes the value n at time t. We 
consider models in which the evolution of these probabilities over time is given by 

f = H(t)p(t) , (4) 

where H(t) is a time-dependent matrix such that at any time t its off-diagonal elements are positive 
and its column sums are zero. This defines a time-inhomogeneous continuous-time Markov chain. 
For some special cases, analytic solutions can be obtained. But in general, where the state-space 
of the Markov chain is finite, numerical algorithms exist that calculate p(t) by making use of 
expansions such as 

p(t + 5t) = H(t)p(t)St + 0(6t 2 ) , (5) 

and accumulating a sufficient number of 8t steps to reach p(t) from p(0). Methods based on this 
direct integration will therefore calculate p(t) in 0(t) operations. 

2.3 The method of Wei and Norman 

In this section, we recall the methodology of Wei and Norman [Tl] , which is applicable to equations 
of the form ^ above. The first step is to look for a decomposition of H(t) 

m 

# (i) = X> ' ( 6 ) 

i=l 

where the Hi are linearly independent matrices obeying 

[Hi,Hj] = H t Hj - HjH l = ^2 ii] k Hk , (7) 

k 

for (in our case real- valued) scalars £jj . Given such matrices, we can form a vector space V = 
spanjiJ;}™ 1 C GL(n, R) such that for a Lie bracket as defined in we have [X, Y] E V for all 
X, Y £ V. This closure under the action of the Lie bracket can be used to look for solutions of the 
form 

p(t) = e 9l{t)Hl ■ ■ ■ e g ™W H ™p(0) =: U{t)p(0) , (8) 

where matrix exponentiation is defined through the power series in the standard way. Using the 
ad operator as defined in ([3]), the Baker-Campbell-Hausdorff formula is 

e x Ye- x = e^Y , (9) 

which will enable us to derive the solution form advertised. Substituting © and ([S]) into (dJ) then 
gives 

J = 5>(Wr(t) P (o) 

(i-l \ (m \ ( 10 ) 

= £*(*) (lie^)^ [He^ P (0). 

*=1 \i=l / \j=i J 

Since this expression holds for any p(0), we can equate the operators acting on p(0), postmultiply 
by U~ x , and repeatedly apply © to obtain 

m ml i—1 \ 

goiW^ = Y,9i(t) n« 93(f)(adflj) H, . (11) 
«=1 i=l \j=l J 
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The precise solution to this equation will depend on the constants £ in ([7]) , however since the Hi are 
chosen to be linearly independent, terms in in front of the same basis matrix can be equated, 
leading to a set of ODEs for the gi(t). 

The usefulness of this method therefore depends on whether appropriate Hi matrices can be 
defined, so that the equations that must be solved for gi(t) are not excessively complex. But 
in the event that gi(t) can be calculated in O(l) rather than, say, 0(t) - which will typically 
be the case if an analytic result is obtained - then the computation of pit) can be achieved in 
0(1) rather than 0(t) through the numerical calculation of the matrix exponentials in (JHJ. Such 
enhanced computational tractability of stochastic models clearly has benefits in the application of 
probability theory to statistical inference, where the speed of evaluation of likelihoods is highly 
important. 



3 Examples 

The primary difficulty in applying the method above to population models is finding an appropriate 
expansion of the form ([5]), since the systems involved are not obviously symmetric. We now turn 
to three examples where an appropriate expansion can be found. In two cases, special initial 
conditions give analytic results that can be checked against other methods to confirm the soundness 
of the approach; and in the third, a significant numerical benefit is observed compared to direct 
integration. 



3.1 A birth-death process 

Suppose we have a time-inhomogeneous birth-death process, characterised by a stochastic variable 
N(t) > taking integer values n, and the transition rates 

n — > n + 1 at rate bit) , 

(12) 

n — > n — 1 at rate nd(t) . 

A biological interpretation of this process would be the survival of juvenile animals, introduced to 
an inhospitable region by seasonal breeding happening at another site, and dying at a rate that 
depends on the climate. Defining components of a vector p n (t) — Pr(7V(i) = n), the Kolmogorov 
equation for this process is 

^ = (b(t)(R-t) + d(t)(L-M)) P . (13) 

The matrices involved are countably infinite in dimension, and are defined implicitly by (1121) 
and (TIB")) . It is also possible to write explicit definitions in terms of the Kronecker delta: 

(l)n,fc = &n,k ; (R)n,k = $n,k+l \ (L) n ,k = (k — l)8 n ^-l 5 (M) n ,k = k5 n ,k ■ (14) 

Clearly, the identity matrix commutes with everything (i.e. [1, X] = for any X) while the other 
matrices obey 

[L, R] = 1 , [A/, R] = R , [L, M] = L . (15) 
We then look for solutions of the form 

p(t) = e ftWie»WH e »W£e*W"j)(0) , (16) 
noting that e 91 ^ 1 = e 9l ^\ Making use of the result (fTTj) together with the algebra ([T5|) gives 



9l (t) = -e v ^ / b(u)e v ^du 
Jo 

92(t)=e v ^ f b(u)e v Mdu , 



9s(t)=e^-l, (17) 
g 4 {t) = -V{t) , 



where T)(t) := / d(u)du 
Jo 
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This provides a solution to the original model, but one that is much simpler if we assume the initial 
condition N(0) = 0, in which case 



Pn{t) = e 9l( 



(18) 



for gi, gi as in (|17p . In this way an infinite-dimensional time-inhomogeneous Markov chain is 
reduced to carrying out the two integrals in (JTTJ) . It is worth comparing this to the 'textbook' 
method for dealing with time-inhomogeneous Markov chains, which is to derive an expression for 
the probability generating function (PGF) [3]. This is done by writing down the Kolmogorov equa- 
tion (|T3|) in component form and substituting in the definition of the PGF, G(s, t) := J2 n sn Pn(t)- 
This gives a PDE for the PGF of 

' (b(t)s-d(t))(s-l)^ , (19) 



dt 

which, given the initial condition A^(0) — 0, gives 



ds 



G(s,i)=exp (s-l)e- v W / b{u)e v{u) du 
l Jo 

The equation (f20|) yields the same solution as (fl~8f above through the standard relation 

1 d n G 



Pn = -f 



! ds r ' 



(20) 



(21) 



s=0 



The effort in deriving the solution (fT8"|) through the PGF and Lie algebraic methods is therefore 
roughly similar; however the intermediate results obtained in each method are likely to be useful in 
different contexts. For example, the PGF in (|2T)]) is likely to be the easiest way to derive moments 
of the process; while the matrix exponential form (flB]) is likely to be useful if the derivative of the 
solution with respect to a parameter of the model is required [15] . 



3.2 Epidemic surveillance 

Consider the following situation. An epidemic is in progress in a population, such that individuals 
are either susceptible to infection, infectious, or recovered and immune. Surveillance of the epi- 
demic is carried out by recruitment of individuals at random from the general population (or more 
realistically through recruitment of individuals in contact with the healthcare system due to non- 
infectious illness) who are tested and determined to be either susceptible, infectious or recovered. 
The epidemic is characterised by a force of infection X(t), which is the rate at which susceptible 
individuals become infectious and for which there are various parametric forms 7], and also by a 
recovery rate 7(i), which is the rate at which infectious individuals recover. A plausible explicit 
choice for these functions is to hold ^(t) constant, and to take X(t) = \oe rt , representing the early 
exponential growth phase that is common to many different epidemics. 

As other authors have found, manipulation of more complex Markov chains is simplified by the 
use of Dirac notation [HJ 03 [5J [T3] . In this formalism, the probability vector for the model described 
above is written 

\p(t))=J2^(S,I\N,t)\S,I) , (22) 
s,i 

where Pi(S, I\N,t) is the probability that from a cohort of size N a time t after the start of the 
epidemic a number S of the cohort is susceptible and a number / is infectious (leaving N — S — I 
recovered individuals). \S,I) is a basis vector, linearly independent of any other basis vector 
with different susceptible and infectious counts. Operators, marked out using a hat O, act on 
basis vectors to give linear combinations of basis vectors. For this system, we need the following 
operators: 



S\S,I) 


= S\S,I) , 


i\s,i) 


= I\S,I) , 


A|5,J) 


= S\S-l,I) , 


p\s,i) 


= I\S,I-1) , 


f\S,I) 


= 5|5-1,J + 1> 
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The action of these operators can be described computationally as follows: S returns the number of 
susceptibles; / returns the number of infectives; A returns the number of susceptibles and depletes 
these by one; p returns the number of infectives and depletes these by one; and f returns the 
number of susceptibles, depletes the susceptible population by one, and increases the infectious 
population by one. The dynamical model is then 

j t IK*)> = {lit) {P-I)+ Kt) (r - S)) W)) . (24) 

Note that while (l2~4"|) does not make use of the operator A, it is necessary to include this to have 
an algebra that is closed under the action of the Lie bracket. The full set of Lie brackets is shown 
in Table [H and the action of the exponential adjoint endomorphism is shown in Table [5J We then 
look for a solution of the form 

\p(t)) =e 9i( t )A e9 2(t)f e9 3WS e94 (t)P e 95(t)/| 2 ,(0)) . (25) 

Going through the same procedure as before gives solution 

gi (t) = e A « - e- A « - j* \(u)e- A ^e r ^~ r ^du^j , 

92(t) = e A « /* \(u)e~ A ^e r ^- r ^du , 
Jo 

g 3 (t) = -A(t), g4t) = l-e- r ^, g 5 (t) = -T(t) , 

for A(t) := / X(u)du , and T(t) := / ^(u)du . 
Jo Jo 

To check this result, assuming \p(0)) — \N,0) and substituting (j2l))) into (f2"5"j) gives 
Pr(SJ\N,t)= { sw{N m s _ 1}l )(*i) S fa) I (l-ni-ir a ) 1 



(26) 



■.N-S-I 



(27) 



where 

TTi = e" A W , n 2 - /* A( H )e- A W e r W-rW rfu . ( 28 ) 
Jo 

This is the form we would expect for this solution; considering each individual's probability of 
remaining susceptible to be tt\ and being infections to be 7T2, these should obey 



dt I TT2 J \ A(t) - 7 (t) M 7T 2 



(29) 



which has solution (f28|) . and the independence of each individual leads to the multinomial distri- 
bution (|27|) for the cohort as a whole. As for the birth-death process, the analytic result obtained 
through Lie algebraic methods can be obtained otherwise, and which method is preferable will 
depend on which further calculations one wishes to undertake. 



3.3 A pure birth process 

Now suppose we have a pure birth process, characterised by a stochastic variable N(t) and transi- 
tion rate 

n — > n + 1 at rate a(t) + nb(t) . (30) 

Special forms of a(t), b(t) have been used to model the formation of social contacts [1 . Defining 
components of a vector p n — Pr(N(t) — n), and assuming that there is a maximum count of 
interest m, such that we only keep track of Pr(iV(t) > m) as the final component of p, we write 
the Kolmogorov equation for this process as 

^ = (a(t)Pi + b(t)Q)p . (31) 
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The generating matrices Pi and Q are implicitly denned by (j3"U|) . however an additional set of 
matrices are required to produce a vector space that is closed under the Lie bracket. These 
matrices take a similar form to P\ and are indexed by integer i. The definition of the matrices 
used, in terms of the Kronecker delta, is 

(Q)n,k = k(S n ,k+l — Sn,k) , (Pi)n,k = <5n,fc+i _ <5n,fc+i-l ■ (32) 

These matrices satisfy the commutation relations 

[P,P,-] = , [Pi, Q] = -iP l+1 + (i - 1)P . (33) 

Note that it is here that the assumption of a finite state space (i.e. the introduction of a maximum 
count of interest above) allows the solution method to work, since otherwise a countably infinite 
number of matrices would be needed. We then look for a solution of the form 

p(t) = e^ p > ■ ■ ■ e /™(*)^ e 9(*)O p (0) , (34) 

and make use of to give the following system of equations: 

A(i) = a{t) , g(t) = b(t) , j t>1 {t) = (i - l)g(t) [/^(t) - £(*)] . (35) 

These equations have solution 

fl(t)=A(t) = f a{u)du , 
Jo 

g{t) = B{t) = [ b(u)du , (36) 
Jo 

/ i>x (t) = e-G- 1 ) 8 **) J* a(u) (e B (*) - e^)^ du . 
Figure [T] shows numerical results for this system for the simple choice 

a(t) = 1 , b(t) = —— , m = 100 , N(0) = . (37) 
Note that for this initial condition, and due to the relations (I33p , we can write the solution as 

Pn{t) = f e ^mm) , (38 ) 

V / n,0 

the evaluation of which which can be seen in Figure[T]to give a significant computational advantage, 
as implemented in EXPOKIT [TT], compared to direct integration of (f3~Tj) through Runge-Kutta, 
as implemented in the MATLAB function ode45, at large times. Perhaps unexpectedly, this is 
seen despite the relatively large value of m. These plots show that, while the ODE solver uses a 
more sophisticated relationship than ([5]) to obtain better than 0(t) performance, it is still much 
more sensitive to model time t than the matrix exponential method. 

There will, of course, be more complex systems where numerical integration via Runge-Kutta is 
impractical, but analytic integration is simple (e.g. if either a(t) or b(t) is a square wave rapidly os- 
cillating between and 1); for these systems the matrix exponential method will further outperform 
Runge-Kutta. But there will also be systems where direct numerical integration is straightforward 
and there is no simply obtained form for fi(t) and g{t), meaning that the matrix exponential 
solution is not useful. 

4 Discussion 

This paper has considered the solution of time-inhomogeneous Markov chains in population mod- 
elling through the use of matrix exponentials. This is done using the method of Lie algebras 
originally developed for applications in physical sciences [14] . In contrast to physical applications, 
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population models are often insufficiently symmetric to write down a well studied Lie algebra. In 
even the relatively simple pure birth process considered, a large number of basis matrices were 
needed to derive a matrix exponential solution; but despite this the exponential solution is useful 
if a derivative with respect to a model parameter is required |15) . and perhaps more importantly 
often has a significant numerical advantage over direct integration of the ODE system [8], Given 
the popularity of computationally intensive inference in modern population models [3], any such 
improvement in numerical efficiency of likelihood evaluation can have important practical benefits. 
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Figure 1: Numerical results for a time-inhomogeneous pure birth process. Top: probability dis- 
tributions at different times. Matrix exponential solutions are shown with markers and direct 
integration via Runge-Kutta is shown with continuous black lines — clearly, these are numerically 
indistinguishable. Bottom left: values of the quantities /j over time. Bottom right: CPU time 
needed to run each method as a function of model time. 
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Table 1: Values of [X, y] for the epidemic model. 
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Table 2: Values of e x ( adX )Y for the epidemic model, for scalar x. 
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